Spatial transcriptome-guided multi-scale framework connects P. aeruginosa metabolic states to oxidative stress biofilm microenvironment

With the generation of spatially resolved transcriptomics of microbial biofilms, computational tools can be used to integrate this data to elucidate the multi-scale mechanisms controlling heterogeneous biofilm metabolism. This work presents a Multi-scale model of Metabolism In Cellular Systems (MiMICS) which is a computational framework that couples a genome-scale metabolic network reconstruction (GENRE) with Hybrid Automata Library (HAL), an existing agent-based model and reaction-diffusion model platform. A key feature of MiMICS is the ability to incorporate multiple -omics-guided metabolic models, which can represent unique metabolic states that yield different metabolic parameter values passed to the extracellular models. We used MiMICS to simulate Pseudomonas aeruginosa regulation of denitrification and oxidative stress metabolism in hypoxic and nitric oxide (NO) biofilm microenvironments. Integration of P. aeruginosa PA14 biofilm spatial transcriptomic data into a P. aeruginosa PA14 GENRE generated four PA14 metabolic model states that were input into MiMICS. Characteristic of aerobic, denitrification, and oxidative stress metabolism, the four metabolic model states predicted different oxygen, nitrate, and NO exchange fluxes that were passed as inputs to update the agent’s local metabolite concentrations in the extracellular reaction-diffusion model. Individual bacterial agents chose a PA14 metabolic model state based on a combination of stochastic rules, and agents sensing local oxygen and NO. Transcriptome-guided MiMICS predictions suggested microscale denitrification and oxidative stress metabolic heterogeneity emerged due to local variability in the NO biofilm microenvironment. MiMICS accurately predicted the biofilm’s spatial relationships between denitrification, oxidative stress, and central carbon metabolism. As simulated cells responded to extracellular NO, MiMICS revealed dynamics of cell populations heterogeneously upregulating reactions in the denitrification pathway, which may function to maintain NO levels within non-toxic ranges. We demonstrated that MiMICS is a valuable computational tool to incorporate multiple -omics-guided metabolic models to mechanistically map heterogeneous microbial metabolic states to the biofilm microenvironment.


Introduction
Biofilms are three-dimensional, dynamic, heterogeneous microbial communities.Emergent biofilm metabolite gradients result in the spatial organization of bacteria with distinct metabolic functions, impacting the fate of the microbial community to expand, cause infection, and resist antibiotics [1,2,3].Gaining a mechanistic understanding of biofilm metabolic spatial organization is limited by current experimental tools, making it difficult to measure the connected intracellular and extracellular metabolic processes across multiple temporal and spatial scales [4].For example, at single-cell resolution, spatial transcriptomic experiments retain the spatial locations of bacterial metabolic states within biofilms [5], but they do not couple that information with high-resolution metabolomics measurements that can reveal mechanistic relationships between metabolic states and metabolite microenvironments [6].To overcome these current experimental limitations, computational tools that bridge mechanisms across temporal and spatial scales can reveal underlying intracellular and extracellular mechanisms and prioritize experiments to test potential therapeutic strategies to treat infectious biofilms [7,8].
One computational approach for simulating intracellular microbial metabolism uses genome-scale metabolic network reconstructions (GENREs).GENREs represent the complex, interconnected metabolic reaction network within a cell using a mathematical description of known gene-protein-reaction relationships and the stoichiometry of associated chemical transformations.GENREs can be interrogated with constraint-based flux-balance analysis (FBA) to simulate flux distributions associated with catabolic and anabolic processes [9].Algorithms are being developed to integrate -omics data into a GENRE to constrain the intracellular metabolic solution space and generate biologically-relevant predictions of cellular metabolism in various metabolic environments [10][11][12].
Although -omics data integration algorithms enable improved prediction accuracy of cellular metabolism, outputs from FBA simulations often represent the steady-state metabolism of a given cell type or species and thus do not capture metabolic heterogeneity across space and time.Dynamic FBA can be used to predict temporal changes in biomass and extracellular metabolite concentrations, but lacks consideration of metabolic spatial heterogeneity [13].To reveal metabolic differences in space, previous efforts integrated spatial transcriptomics data from healthy and cancerous tissue regions into a human GENRE, but lacked predictions of dynamic cell-cell and cell-environment interactions to predict disease dynamics [14].Despite these efforts, few computational frameworks have simulated -omics-integrated GENREs in both spatial and temporal dimensions [15], which is important to mechanistically predict pathogenesis and therapeutic outcomes in addition to other biological process.
To simulate emergent spatiotemporal metabolic heterogeneity, computational frameworks have coupled GENREs with an agent-based model (ABM), which simulates individual cell behavior, and a reaction-diffusion model solved with partial differential equations (PDEs) that predicts extracellular metabolite concentrations.Previous 2D multi-scale frameworks, BacArena, MATNET, and COMETS, represent individual agents as a single cell or a population of cells, in which the simulation frameworks predicted emergent metabolic heterogeneity because each agent's GENRE was constrained by local nutrient fluxes within a heterogenous nutrient environment [8,16,17].Despite these efforts, few multi-scale frameworks are easily extendable for cells to dynamically adopt different -omics-integrated GENREs, which may improve predictions of metabolic processes controlled by gene regulation mechanisms.In addition, multiscale frameworks have not incorporated GENREs integrated with spatial transcriptomic data, which can capture metabolic heterogeneity at single-cell spatial resolutions.For example, because the 3D multi-scale framework ACBM implemented a GENRE integrated with population-level transcriptomic data, this framework was not likely to capture potential biofilm metabolic heterogeneity measured at single-cell resolutions [15].
In this work, we present an extendable multi-scale computational framework that couples multiple -omics data-integrated GENREs, an ABM, and metabolite reaction-diffusion PDEs.We refer to this framework as a Multi-scale model of Metabolism In Cellular Systems (MiM-ICS).MiMICS is an open-source Java-and Python-based framework.MiMICS is extendable to simulate in 2D and 3D, and to represent individual agents as a single cell or a population of cells.A key feature of MiMICS is the ability to incorporate multiple -omics data-integrated GENREs, which can represent unique metabolic states.As a result of the corresponding integrated -omics data, each metabolic model state may predict different parameter values that alter the extracellular environment, such as nutrient uptake or toxic byproduct secretion.MiMICS allows for the user to incorporate multiple GENREs integrated with -omics data measured at the single-cell or population-scale level.While MiMICS was designed to integrate spatially resolved transcriptomics data, MiMICS could be used to integrate global transcriptomics data that was measured in various metabolic conditions.MiMICS can execute biologically-relevant ABM rules for cellular agents to choose from the -omics data-integrated GENREs to simulate metabolism.Simple mechanistic rules were used for a cell to switch metabolic model states related to the cell's extracellular metabolic environment, effectively representing the transcription of metabolic genes regulated by the metabolic environment.Future studies could directly couple MiMICS with a gene regulatory network, which predicts gene transcription regulated by the metabolic environment [18].
As an initial biological test case to demonstrate its utility, MiMICS was applied to simulate emergent metabolic heterogeneity within a 3D Pseudomonas aeruginosa biofilm observed by a recent spatial transcriptomic study [5].P. aeruginosa is an opportunistic pathogen that can cause deadly biofilm infections in the lungs of patients with cystic fibrosis and COVID-19 [19,20].The published spatial transcriptomic study revealed microscale spatial organization of aerobic, denitrification, and oxidative stress metabolic states within a P. aeruginosa PA14 biofilm [5].As proposed by Dar and co-workers [5], genes related to denitrification metabolism, an anaerobic respiration process, were hypothesized to be upregulated in anoxic PA14 biofilm regions [5,21,22].Likely secreted by denitrification cells, the cytotoxic denitrification intermediate nitric oxide (NO) was also hypothesized to upregulate oxidative stress genes in nearby PA14 biofilm cells [23].However, the experiment lacked a quantitative and mechanistic mapping of the location of the cell and its metabolic state to the surrounding metabolic microenvironment [5].Thus, herein, an established algorithm (RIPTiDe) [11] was used to integrate the published P. aeruginosa PA14 biofilm spatial transcriptomic dataset into a previously curated P. aeruginosa PA14 GENRE that generated four unique PA14 metabolic model states.The metabolic model states captured differences in aerobic and anaerobic denitrification metabolism, and revealed denitrification subpopulations that secreted the cytotoxic metabolite NO.This latter metabolic model state was crucial to predict a NO secretion rate that was passed to the extracellular reaction-diffusion model in MiMICS to generate a NO biofilm microenvironment that induced oxidative stress.Agents decided which metabolic model state to simulate their intracellular metabolic processes based on a combination of stochastic rules and metabolite sensing rules, the latter considering oxygen and NO levels in the agent's local environment.MiMICS predicted microaerobic and variable NO microenvironments emerged within biofilm regions, resulting in microscale patches where cells heterogeneously used denitrification and oxidative stress metabolism.Due to cells sensing extracellular NO signals, MiMICS revealed the dynamics of cell populations heterogeneously regulating reactions in the denitrification pathway, which may function to maintain NO biofilm concentrations within non-toxic ranges.As demonstrated with this P. aeruginosa biofilm test case, we believe MiMICS is a promising computational tool that can use multiple -omics data-integrated metabolic models, and mechanistically simulate and map heterogeneous microbial metabolic states to the biofilm microenvironment.

Results & discussion
MiMICS is an extendable computational framework executed in Python and Java to simulate metabolism in 2D and 3D microbial communities.MiMICS couples a genome-scale metabolic network reconstruction (GENRE) with the established platform Hybrid Automata Library (HAL) [24], which contains an agent-based model (ABM) and a continuum-scale reaction-diffusion model (Fig 1).To our knowledge this is the first multi-scale metabolic framework to interface with the COBRApy Python package, which is becoming increasingly common to simulate and integrate -omics data into a GENRE [11,25,26].In addition, MiMICS is the first framework to couple an intracellular metabolic model with HAL, and simulate 3D microbial biofilms using HAL, which can be challenging due to the small microbial length scales that can cause instability in PDE solvers [24].MiMICS offers the user the ability to input multiple -omics data-integrated metabolic models, which can represent distinct metabolic states and yield different metabolite uptake or secretion rates that are passed to the extracellular reactiondiffusion model.Individual agents decide which metabolic model state to execute based on mechanistic rules input by the user, such as agents' sensing of their local metabolite concentrations.
In this study, an individual 2 μm x 2 μm x 2 μm agent represented a single-cell P. aeruginosa bacterium, which existed on a 230 μm x 230 μm x 40 μm three-dimensional grid, corresponding to experimental microscopy dimensions [5].The continuum-scale reaction-diffusion model simulated oxygen, nitrate, nitric oxide, and glucose concentrations within an equivalently sized three-dimensional metabolite grid solved with partial-differential equations [24].Agents were randomly initialized at t = 0hrs, and MiMICS simulated biofilm growth for ten hours, replicating the experimental biofilm growth period [5].At each five-minute simulation time step, each model component of MiMICS was executed to update agent properties and metabolite concentrations (Fig 1).For each agent, the agent's biomass and metabolite concentrations from the continuum-scale grid corresponding to the agent's location were converted to metabolite uptake fluxes used to constrain the agent's intracellular metabolic model (Fig 1).Constraint-based flux-balance analysis was used to optimize each agent's metabolic model to predict a biomass growth rate, as well as metabolite secretion and uptake rates.The biomass growth rate was passed to the ABM to update an agent's biomass (Fig 1).Each agent's predicted metabolite secretion and uptake fluxes were passed to the reaction-diffusion model to update the metabolite concentrations at each agent's location (Fig 1).In the ABM, bacteria agents performed cell division, moved via motility, and performed cell mechanical behaviors to prevent cell overlap (Fig 1).Simulation outputs such as agent locations, agent intracellular metabolic fluxes and metabolite concentrations were reported at desired simulation time steps.

Construction of spatial transcriptome-guided PA14 metabolic models
A key feature of MiMICS allows users to input multiple -omics-guided metabolic models, which may differ in predicted intracellular and extracellular metabolic parameter values, into a multi-scale metabolic framework (Fig 1).A set of four transcriptome-guided PA14 metabolic models was constructed by integrating an established PA14 GENRE [12] with a published spatial transcriptomic P. aeruginosa PA14 biofilm dataset [5] using the established algorithm RIP-TiDe [11].In addition, the PA14 GENRE was constrained on synthetic cystic fibrosis sputum medium (SCFM), replicating SCFM used in the experiment [5].To avoid the computationally expensive generation of transcriptome-guided metabolic models for all ~292,000 PA14 biofilm cells measured, transcriptomic data representative of four unique metabolic states was extracted using the UMAP Leiden clustering method [5] and manual selection (refer to Methods).This approach generated four unique transcriptome-guided PA14 metabolic models.Each represented a unique metabolic state within the ten-hour PA14 biofilm, designated as (1) aerobic respiration, (2) denitrification, (3) denitrification + nitric oxide (NO) secretion, and (4) oxidative stress.
The aerobic respiration metabolic model, integrated with transcriptomic data of cells with high expression levels of ccoN1, which encodes the primary aerobic respiration oxidase [27], predicted high levels of oxygen uptake flux (Fig 2D ).A denitrification PA14 metabolic model was generated from cells with high expression levels of genes encoding denitrification reductases, narG, napA, nirS, norB, and nosZ (Fig 2A), which convert nitrate to nitrogen to perform respiration in low oxygen environments [22] (complete PA14 GENRE denitrification pathway shown in Fig 2C).The denitrification metabolic model accurately predicted flux through the denitrification reductase reactions, resulting in characteristic nitrate uptake and nitrogen secretion fluxes (Fig 2B and 2D) [28].A denitrification +NO secretion PA14 metabolic model was generated from cells with high narG and nirS expression, respectively encoding nitrate and nitrite reductase, but low norB expression, encoding NO reductase.This expression profile, considered to be limited by the norB expression, resulted in a metabolic model state with nitrate uptake and NO secretion flux (Fig 2D).When simulated by individual agents in MiM-ICS, this NO secretion rate was passed from the agent's intracellular scale (GENRE) to extracellular scale (reaction-diffusion model), which is an essential step to produce the NO biofilm microenvironment that induces oxidative stress.Demonstrating the impact of transcriptomic data integration, the PA14 metabolic model without transcriptomic data integration, termed transcriptome-free in Fig 2D, predicted oxygen and nitrate uptake, but did not predict NO secretion denitrification metabolic processes (Fig 2D).
Lastly, an oxidative stress PA14 metabolic model was generated from cells with high expression of the gene katA, which encodes the antioxidant enzyme Catalase A [23] (Fig 2A).Accordingly, the PA14 GENRE was updated with NO cytotoxic and antioxidant Catalase A intracellular reactions [23].Briefly, to mitigate cytotoxic NO binding to intracellular iron, Catalase A binds to NO, acting as a NO sink supplementary to NO degradation by NorB [23] (Fig 2C).Correspondingly, the oxidative stress metabolic model predicted NO was consumed (Fig 2D ) and preferentially degraded by NorB for biomass synthesis before being diverted through the antioxidant Catalase A reaction to mitigate NO toxicity [23] Constrained on replete SCFM, the aerobic PA14 metabolic model state predicted a biomass growth rate similar to an experimental growth rate of aqueous P. aeruginosa in SCFM [29] (S2 Fig) .The denitrification +/-NO secretion PA14 metabolic models predicted the lowest biomass growth rates indicative of low oxygen conditions [30] (S2 Fig) .In addition to the predicted NO secretion rates, the predicted oxygen, nitrate, and glucose exchange flux rates (Fig 2D ) were provided as inputs to update the corresponding metabolite concentrations in the reaction-diffusion model.Due to the differences observed in gltB expression, which encodes a glucose binding protein, among metabolic model states (Figs 2A and Fig 2B), and previous reports relating denitrification metabolism and glucose uptake [31,32], glucose was chosen as the carbon source to simulate in the reaction-diffusion model.

Mechanistic incorporation of transcriptome-guided PA14 metabolic models into MiMICS
MiMICS employs ABM rules for an agent to select a metabolic model state from a set of transcriptome-guided metabolic models input by the user.In this study, a combination of mechanistic ABM rules informed from literature and the spatial transcriptomic dataset were implemented into MiMICS (Table 1).P. aeruginosa has been observed to upregulate genes related to aerobic and denitrification metabolism in aerobic and deplete oxygen conditions, respectively [22].In addition, in the presence of extracellular NO, P. aeruginosa increases expression of katA, encoding the antioxidant Catalase A [23].Thus, to decide between the four PA14 metabolic model states, agents compared their local oxygen and NO concentrations to respective concentration thresholds, [O 2 ] t and [NO] t (Table 1).The parameter value for [NO] t was obtained from literature as the extracellular NO concentration that induced katA expression [23].The parameter value for [O 2 ] t was fit to experimental outputs (S4 Fig) .An oxic O 2 threshold, 0.21 mM, had the smallest MiMICS model error, which suggests complete oxygen depletion was not essential to induce denitrification metabolism.Indeed, previous studies observed P. aeruginosa used denitrification in microaerobic conditions (~0.05 mM oxygen) as a possible supplementary or competitive respiration strategy to aerobic respiration [28].The parameterized [O 2 ] t value in the oxic range can be expected because the height (~10 μm) of the ten-hour biofilm was not large enough to predict significant oxygen depletion [33].In the case that lower oxygen concentrations were present in the biofilm, other mechanisms that deplete oxygen near the biofilm, such as an oxygen boundary layer [34] or oxygen consumption by the planktonic phase above the biofilm may be present.These mechanisms serve as updates to simulate metabolite transport in upcoming MiMICS versions.Altogether, MiMICS suggests an oxygen transition point emerged within the PA14 biofilm that induced denitrification metabolism, but oxygen depletion mechanisms in regions outside of the biofilm may have been present.
A stochastic parameter, R n , generated only for agents in the denitrification-inducing low oxygen and low NO condition, was used for agents to select between the two denitrification metabolic model states, one with and one without predicted NO secretion flux (Table 1).The R n threshold was informed from experimental data, which suggested 6% of denitrification cells exhibited the norB-limiting gene expression profile associated with NO secretion (S2 Fig).This R n parameter suggests cells stochastically expressed nirS or norB, which encode nitrite reductase and NO reductase respectively.Indeed, stochastic expression of denitrification genes has been observed in other bacterial species [35].In addition, experimental observations of P. aeruginosa biofilms after ten hours showed high expression levels of pilA (S1 Fig), which encodes type IV pili protein (PilA) that facilitates surface motility and shapes biofilm structure [5,36].As pili synthesis reactions were not in the current PA14 GENRE [37], simple surface motility ABM rules were incorporated to recapitulate the PA14 biofilm structure and total cell count (S3 Fig).
Table 1.Mechanistic ABM rules for P. aeruginosa agents to choose a PA14 metabolic model state.P. aeruginosa agents compared their local oxygen and nitric oxide concentration to respective metabolite thresholds.A stochastic parameter, R n , determined the agent's decision between the denitrification +/-NO secretion metabolic models.

PA14 metabolic models
Oxygen

Transcriptome-guided MiMICS predicted microscale metabolic heterogeneity and NO microenvironment in PA14 biofilm
MiMICS was simulated with the mechanistic ABM rules controlling agent execution of a transcriptome-guided PA14 metabolic model state (Table 1) (referred to as transcriptome-guided MiMICS simulation).Qualitatively, in comparison to the experiment, the transcriptomeguided MiMICS simulation accurately predicted microscale, spatially confined biofilm niches in which cells heterogeneously upregulated denitrification and oxidative stress metabolic processes (Fig 3).In the simulation, these niches were located near the center of the biofilm colony, where microaerobic and variable NO concentrations were predicted (Fig 3).Thus, MiMICS simulations suggested that variable NO signal concentrations within microaerobic biofilm regions resulted in a heterogeneous population of cells using denitrification and oxidative stress metabolism co-existing within the same microscale niche.Qualitatively, experimental PA14 biofilms observed greater spatial dispersion of cells expressing denitrification and oxidative stress genes (Figs 3 and S5), suggesting other mechanisms may regulate denitrification gene expression, such as stochastic expression of the denitrification gene narG encoding nitrate reductase [35].
To demonstrate the advantages of incorporating multiple transcriptome-guided metabolic models, MiMICS was simulated with the PA14 GENRE unconstrained by transcriptomic data (referred to as transcriptome-free MiMICS), which is the standard practice for current multiscale metabolic frameworks [8,16,17].Transcriptome-free MiMICS inaccurately predicted a homogeneous biofilm population with active flux through all intracellular denitrification metabolic reactions (Figs 3 and S5).In addition, transcriptome-free MiMICS did not predict extracellular NO in the biofilm microenvironment, resulting in agents lacking flux through the NO-induced oxidative stress reaction encoded by katA (Figs 3 and S5).

Transcriptome-guided MiMICS improved prediction accuracy of heterogenous PA14 biofilm metabolism
Quantitative validation of MiMICS was first performed by comparing the experimental percentage of total cells expressing a gene to the simulation's percentage of total agents with active flux through metabolic reactions encoded by that respective gene (Fig 4).Compared to transcriptome-free MiMICS that predicted 100% of cells with active denitrification reaction flux, transcriptome-guided MiMICS accurately predicted reduced cell populations with active flux through denitrification narG-, napA-, nirS-, norB-, and nosZ-encoded reactions (Fig 4).
As a result of the ABM rules enforced in transcriptome-guided MiMICS (Table 1), a heterogenous agent distribution simulating the denitrification and oxidative stress metabolic model state emerged within the spatially variable NO biofilm microenvironment ( Interestingly, transcriptome-guided MiMICS did not accurately predict larger cell populations expressing nirS (20.4% of cells in experiment, 10% in simulation) relative to cell populations expressing narG, norB, and nosZ (specific percentages reported above) (Fig 4).While NO has been observed to also upregulate nirS expression 32 , hypothesized to bind NO required for NO transport to NorB for degradation [38], the current PA14 GENRE does not include this  Potential inaccuracies in either the -omics data-integration method or the MiMICS model components rendered the transcriptome-guided MiMICS incapable of predicting the relatively low cell proportions utilizing gshB and aceA, which encode glutathione synthetase and isocitrate lyase, respectively, or cysA, which encodes a sulfate transport protein (Fig 4).In addition, neither MiMICS transcriptome version captured heterogeneous expression of genes related to virulence factor synthesis, which was expected as virulence factor synthesis likely does not contribute to the biomass synthesis GENRE objective function used during transcriptome data integration [37].Overall, compared to the standard practice of using a transcriptome-free multi-scale framework, transcriptome-guided MiMICS improved prediction accuracy of the heterogeneous cell distributions with upregulated denitrification and oxidative stress metabolism.
To investigate the impact of the ability for cellular agents to dynamically adopt different metabolic model states, a MiMICS scenario was run with cells that remained fixed in a randomly initialized metabolic state, called 'Fixed metabolic state MiMICS' (S9 Fig) .Fixed metabolic state MiMICS predicted homogeneous metabolic state niches localized to the initial location of the respective metabolic state and limited spatial mixing of different metabolic states.In contrast, the dynamic switching ability in MiMICS promoted microscale niches of cells with a heterogeneous distribution of metabolic states localized to the biofilm center.This result confirms that dynamic metabolic state adaptation to the metabolic microenvironment is essential to recapitulate the confined microscale niches of heterogenous metabolic states observed in the experimental biofilm.

MiMICS captured spatial relationships of intracellular metabolism in PA14 biofilm
To quantitatively evaluate the validity of transcriptome-guided MiMICS predictions in space, a bulk neighborhood spatial correlation analysis between genes was performed, similar to the analysis of experimental data by Dar and co-workers [5]  gene pairs.Specifically, MiMICS accurately predicted that denitrification genes (i.e.narG, nirS, norB, nosZ) were positively correlated with one another, and positively correlated with the oxidative stress gene katA (Fig 5A and 5B).Conversely, MiMICS predicted napA, which encodes nitrate reductase, was positively correlated with the remaining denitrification genes narG, nirS, norB, and nosZ.However, this correlation was not observed in the experiment.Upon closer inspection, MiMICS predicted patches in the biofilm where neighboring cells utilized napA-, nirS-, norB-, and nosZ-encoded reactions (S6 Fig) .In contrast, cells expressing napA in the experiment were more dispersed compared to more spatially confined biofilm regions of cells expressing nirS, norB, and nosZ (S6 Fig) .MiMICS prediction discrepancies in the napA spatial correlations suggest alternate mechanisms exclusively modulate napA expression, such as extracellular phenazine secretion [39].Indeed, phenazine synthesis genes phzE1 and phzM were expressed in the ten-hour PA14 biofilm (Fig 4).In addition, MiMICS accurately predicted denitrification and oxidative stress genes (i.e.narG, nirS, norB, nosZ, katA) were anticorrelated with the carbon metabolism genes gltB and sucC (Fig 5A and 5B).One aspect MiMICS did not accurately predict was the spatial correlation of atpA, encoding ATP synthase, with all other genes (Fig 5B ), which motivates improvements in future transcriptome-guided MiMICS simulations.

MiMICS revealed functional dynamics of denitrification and oxidative stress metabolism in PA14 biofilm
This spatial transcriptomics experiment required biofilm fixation, preventing in-situ dynamic monitoring of gene expression within the biofilm sample.Thus, this spatial transcriptomic method was unable to characterize the temporal events that underline the PA14 biofilm gene expression distributions and spatial patterns.In contrast, transcriptome-guided MiMICS dynamic simulations were used to monitor and quantify spatiotemporal shifts of PA14 biofilm metabolism connected to the metabolite microenvironment.
Agents were randomly initialized at t = 0hrs, nucleating the simulated PA14 biofilm which was grown over a period of ten hours (Fig 6).After approximately nine hours of simulated PA14 biofilm growth, agents exposed to the emergent microaerobic biofilm environment (i.e.local oxygen below the [O 2 ] t value) switched from the aerobic to denitrification metabolic model state (Fig 6A and 6B).In accordance with the aerobic and denitrification metabolic model states (Fig 2A ), cell populations with active flux through narG-, nirS-, norB-, and nosZencoded denitrification reactions increased over time (Fig 6A and 6B).The stochastically-chosen agents that executed the norB-limiting denitrification metabolic model (Table 1) secreted NO, generating the NO biofilm microenvironment (Fig 6A).The maximum level of NO, ~4 μM, below toxic levels [23], was localized in the biofilm center (Fig 6A).At ten hours, agents that sensed local extracellular NO above the [NO] t value switched to the oxidative stress Due to the ABM rules that were enforced (Table 1), a heterogeneous distribution of cells using the denitrification and oxidative stress metabolic models emerged within microaerobic and spatially variable NO biofilm microenvironments.Only the denitrification +/-NO metabolic model states predicted active flux through narG-and nirS-encoded intracellular reactions, which are upstream of NO synthesis and required for NO secretion (Figs 2B and 6A).Conversely, denitrification and oxidative stress metabolic model states predicted active flux through norB-and nosZ-encoded intracellular reactions, which are downstream of NO synthesis and used for NO degradation (Figs 2B and 6A).As a result of these intracellular differences among metabolic model states, fewer cells in the biofilm were predicted to use narG-and nirSencoded reactions compared to cells using norB-and nosZ-reactions (Fig 6B ).
Relating these MiMICS dynamic predictions to possible biological functions, previous studies similarly observed cells downregulated narG expression and upregulated expression of nirS and norB, the latter encoding a NO degradation reductase, in presence of NO [40,41].This NO-induced regulation of different genes in the denitrification pathway was hypothesized to arrest NO synthesis, promote NO degradation, and maintain extracellular NO within nontoxic concentrations [40,41].Functioning similarly, MiMICS predicted NO signaled agents to reduce flux through narG and nirS-encoded NO synthesis reactions, which decreased NO synthesis and secretion (Fig 6).Due the decreased NO secretion, and extracellular NO gradients dissipated at 9.9 hours in the biofilm colony shown in Fig 6A .In contrast, possibly to promote degradation of extracellular NO, agent populations with active flux through norB-and nosZencoded NO degradation reactions continued to increase (Fig 6).
Next, possibly due to the reduced NO concentrations and expanding microaerobic biofilm regions, agents switched to a denitrification +/-NO secretion metabolic model state, again promoting a NO-rich biofilm microenvironment at 10 hours in the biofilm colony in Fig 6A .Thus, MiMICS predicted an oscillating, non-toxic NO biofilm microenvironment (Fig 6).This prediction is supported by previous studies that have measured temporal oscillations in extracellular NO in low oxygen P. aeruginosa cultures performing denitrification [42,43].Furthermore, temporal oscillations in expression of mRNA encoding denitrification reductases has been observed in E. coli [44].In MiMICS, the predicted oscillations in extracellular NO were a result of cellular agents differentially regulating intracellular reactions in the denitrification pathway in response to extracellular oxygen and NO, which may function to maintain a nontoxic NO biofilm microenvironment.Future simulations can systematically perturb genes and parameters in MiMICS to test their effect on predicted outcomes, such as NO metabolism in the biofilm.For example, in silico gene knock-outs, such as ΔnirS that increases NO secretion, or varied extracellular nutrient concentrations may be screened for increased NO-induced P. aeruginosa cell death as a potential therapeutic strategy.Simulation readouts, such as biofilm structures, gene expression, or number of live/dead cells, can be quantified and compared to the experiment.

MiMICS computational performance
Performing GENRE simulations for individual cellular agents is advantageous to predict emergent metabolic heterogeneity within a multicellular system but can be computationally expensive.Execution of one MiMICS simulation time step on one central processing unit (CPU) resulted in a runtime on the order of 10 minutes for 10,000 cellular agents.To improve this MiMICS runtime, GENRE simulations for each individual agent were divided across multiple CPUs using parallel computing, implemented with Python Multiprocessing.This parallel computing strategy reduced MiMICS runtime by an order of magnitude for 10,000 cellular agents (S7 Fig).This computational performance is similar to simulation runtimes reported by the BacArena framework [17] for 1000 agents, where MiMICS is the first framework to report efficient runtimes for up to 10,000 single-cell agents.The number of metabolic model states input into MiMICS did not considerably impact MiMICS runtime (S7 Fig).

Conclusion
In this work, a novel multi-scale computational framework, MiMICS, was developed and presented.MiMICS is an open-source computational framework that interfaces Python-based metabolic models with the established Java-based Hybrid Automata Library agent-based model and reaction-diffusion model platform.MiMICS is extendable for users to input multiple -omics data-integrated metabolic model states that yield different metabolic parameter values that are then passed to the extracellular models.
Used as a test case to demonstrate the framework's utility, MiMICS simulated the connected, multi-scale metabolic processes controlling microscale metabolic heterogeneity observed in a P. aeruginosa PA14 biofilm.MiMICS applied mechanistic ABM rules for P. aeruginosa agents to choose a spatial transcriptomic-guided PA14 metabolic model state.As a result, MiMICS predicted microscale niches where cells heterogeneously upregulated denitrification and oxidative stress metabolism mapped to a microaerobic and variable NO microenvironment.Transcriptome-guided MiMICS revealed the spatial dynamics of cells heterogeneously regulating reactions in the denitrification pathway, possibly to maintain a non-toxic NO biofilm microenvironment.Future MiMICS perturbations could seek to promote NO-induced cell death to inform potential PA14 biofilm therapeutic strategies.While not supported by the experimental dataset that was measured at widely spaced time points, NO-induced biofilm dispersal [45] could be added to the MiMICS mechanistic rule set to simulate dynamic biofilm restructuring.
In this work, the single-cell spatial transcriptomics dataset was used to calibrate the metabolic model component of MiMICS, and MiMICS predictions were tested with metrics at the multicellular level (e.g.gene neighborhood Pearson correlation).To calibrate the metabolic model with -omics data independent of the test dataset, future MiMICS applications could integrate global or single-cell transcriptomics data that was measured independently in relevant metabolic conditions (e.g.varying oxygen concentrations).Future work could also extend MiMICS to explore other mechanisms, such as quorum-sensing signals, that regulate heterogeneous metabolism in single-or multi-species biofilms.Overall, we believe MiMICS is a valuable computational tool to integrate -omics data and elucidate mechanisms that control metabolic heterogeneity, which can promote virulence and antibiotic tolerance in microbial biofilms.

Genome-scale network reconstruction (GENRE)
Pseudomonas aeruginosa PA14 metabolism was simulated with a published P. aeruginosa PA14 genome-scale metabolic network reconstruction (GENRE) (iPau21) [12,37].All GENRE simulations were performed in Synthetic Cystic Fibrosis Sputum Medium (SCFM).A few reactions were updated in the PA14 GENRE to represent the denitrification, fermentation, nitric oxide secretion, and oxidative stress metabolic processes observed in experiments.PA14 GENRE simulations in anaerobic + nitrate were initially infeasible (S8 Fig), which did not reflect expected denitrification biomass synthesis processes [22].In anaerobic SCFM conditions, ubiquinol-9 (UQ 9 ), a preferred respiration cofactor required in the PA14 GENRE biomass objective function [37], synthesis was found to be blocked which prevented biomass growth (S8 Fig) .Three oxygen-independent hydroxylation reactions of UQ 9 precursor metabolites were added with H 2 O as the reactant in the place of oxygen, a hypothesized oxygen-independent reaction [46,47,48].Oxygen-dependent UQ 9 synthesis reactions were set to be irreversible to prevent the generation of intracellular oxygen.To promote feasible L-arginine fermentation, reactions encoded by arcB and argH were set as reversible and irreversible, respectively, in agreement with the expected arginine deiminase pathway [49].These updates to the PA14 GENRE promoted feasible biomass growth rates and flux through the respective reactions in anaerobic conditions (S8 Fig) .Exchange and transport reactions were added for the denitrification intermediates nitric oxide (NO) and nitrous oxide (N 2 O).NO toxicity was represented as a sink reaction for NO and iron [23] and an oxidative stress KatA reaction was represented as a sink reaction for NO [23].Of note, when cells were exposed to an extracellular NO microenvironment in MiMICS, NO uptake was enforced by fixing the GENRE upper bound for NO exchange at a negative flux value, which was calculated from the cell's local extracellular NO concentration (Eq 1).
In MiMICS, each agent's PA14 metabolic model was optimized for a biomass synthesis objective function using flux-balance analysis (FBA) [9].Each agent's metabolic model exchange reaction bounds were constrained by the agent's local oxygen, nitrate, glucose, and nitric oxide concentrations.The agent's local patch metabolite concentrations and metabolic model lower bound of exchange reaction flux were converted with the following equation: Where f M is the metabolite flux that set the metabolic model exchange flux (mmol/(gdWt *hr), c M is the agent's local patch metabolite concentration (mM), v patch is the patch volume unoccupied by the bacteria (L), dt ex is the metabolite uptake time step (estimated as 0.05s in Table 2), and b is the agent biomass (g).v patch was assumed as 1e -16 L.

Construction of transcriptome-guided PA14 metabolic models
The PA14 GENRE was integrated with a published spatial transcriptomics PA14 biofilm dataset [5] using the established RIPTiDe algorithm [11].The dataset consisted of the mRNA expression of ~292,000 PA14 biofilm cells grown in SCFM, and fixed at either 10hr (n = 7) or 35hr (n = 3) growth time points [5].For data integration into the PA14 GENRE, gene expression values were only used for the 47 of the 105 genes measured in the experiment that also encoded reactions in the PA14 GENRE (S1 Table ).For data integration into the PA14 GENRE, lower bounds for exchange reactions were constrained according to SCFM concentrations (S2 Table ), converted using Eq 1.In addition, upper bounds for exchange reactions were set to +1000 to allow for unconstrained production of exchange metabolites.A UMAP Leiden clustering analysis (scanpy v.1.7.0), data preparation and parameters described elsewhere [5], identified distinct gene expression cellular states within the spatial transcriptomic PA14 biofilm dataset (S1 Fig) .Metabolic states for each cluster were assigned based on the highest ranked genes within the cluster.Median gene expression was extracted for the top 9 clusters, which captured 91% of biofilm cells, that were integrated into the PA14 GENRE using RIPTiDe [11] to generate 9 UMAP-identified PA14 metabolic model states.A fractional growth rate of 0.7 was used in RIPTiDe.Reactions not shared by all UMAP-identified metabolic model states, termed non-consensus reactions, and predicted exchange metabolite fluxes were compared and grouped based on similarity to identify unique 10-hour and 35-five PA14 metabolic state models (S1 Fig) .Due to the wide diversity of metabolic states between the 10hr and 35hr time points, the initial scope of this study was on metabolic states emerging within PA14 biofilms grown for 10 hours.This UMAP-informed method generated the 10hr biofilm aerobic PA14 metabolic model state that was input into MiMICS.
A manual selection method was used to extract gene expression values of denitrification and oxidative stress metabolic subpopulations that were not identified in the UMAP Leiden clustering analysis.The ~22,000 denitrification cells classified in the clustering analysis were categorized into having high expression of narG and/or napA (S2 Fig) .The narG-and napAexpressing denitrification-classified cells were further classified into possessing a non-limiting or limiting denitrification pathway (S2 Fig) .Cells with a non-limiting denitrification pathway did not possess a rate limiting denitrification transcriptional step (i.e.expression of narG < = nirS < = norB < = nosZ).Cells categorized with limited denitrification metabolism had expression of a gene encoding a denitrification reductase that was larger than the expression value of the gene encoding the subsequent denitrification reductase gene.From cells classified with non-limiting and limiting denitrification metabolism, the median gene expression was extracted and integrated into the PA14 GENRE using RIPTiDe.The fraction of the optimal biomass was set as 0.55 used in RIPTiDe for these denitrification subpopulations [30].The method generated four PA14 denitrification metabolic model states: no limiting denitrification step, nirS-limiting, norB-limiting, and nosZ-limiting.Predicted exchange fluxes of each of these denitrification limiting metabolic model states are shown in S2 Fig.This method generated the denitrification +/-NO secretion metabolic model states that were input into MiMICS.
An oxidative stress PA14 metabolic model state was generated by integrating the PA14 GENRE with median gene expression values extracted from the ~22,000 denitrification-classified cells with high katA expression.To simulate the exposure of an NO-induced oxidative stress environment, during transcriptome integration into the PA14 GENRE, the upper bound for the NO exchange reaction was set to a negative NO flux value (converted from 20 μM NO) to enforce NO uptake.The fractional optimal biomass was 0.7 used in RIPTiDe.This method generated the oxidative stress PA14 metabolic model that was input into MiMICS.

Metabolite reaction-diffusion model
Metabolite concentrations were simulated with the built-in HAL partial differential equation (PDE) alternating direction implicit (ADI) method solver.In HAL, the metabolite grids were defined as 3-dimensional 230 μm x 230 μm x 40 μm grids divided into 2 μm x 2 μm x 2 μm patches (Table 2).Metabolite grids were constructed for extracellular oxygen, nitrate, nitric oxide, and glucose.
Metabolite concentrations were initialized in each patch based on a Synthetic Cystic Fibrosis Sputum Medium (SCFM) recipe (S2 Table) [29].The initial condition for the metabolite PDEs was assumed to be spatially uniform: cðx; y; z; t ¼ 0Þ ¼ c initial where c initial was the initial metabolite concentration (mM).
Before metabolite diffusion simulations, a metabolite reaction term was discretely applied at each agent's respective grid location, representing the local metabolite consumption or production by that agent.The reaction term consisted of first converting the agent's metabolic flux, predicted by the agent's metabolic model, to a concentration (Eq 1), which was then passed to update the metabolite concentration at the agent's grid location (Eq 2).
where c f is the updated metabolite concentration (mM) at the agent's location, c 0 is the initial metabolite concentration (mM) at the agent's location, and Δc is the agent's metabolic model predicted metabolic flux converted to a concentration (mM).With a consideration towards the thousands of agents, each with multiple predicted metabolic reaction rates, this discretely applied metabolite reaction rate approach was used to improve computational efficiency.Next, diffusion was simulated using HAL's ADI diffusion equation solver method.The governing equation for metabolite diffusion was as follows: where c is the metabolite concentration (mM), t is time (s), x is the x-dimension, y is y-dimension, z is the z-dimension, and D is the metabolite diffusion coefficient (cm 2 /s).
No flux boundary conditions were applied in the z-dimension: @c @z j z¼0 ¼ @c @z j z¼L ¼ 0 where z = 0 was at the glass slide and z = L was at the air-aqueous interface.Periodic boundary conditions in the x-and y-dimensions were applied.The biofilm region was defined as R biofilm .Regions outside of R biofilm assumed constant metabolite concentrations.This assumption was attributed to higher metabolite diffusion coefficients occurring in the aqueous phase compared to metabolite diffusion coefficients the biofilm.In addition, the aqueous nutrient-rich media was refreshed every four hours in the experiment [5], so constant nutrient concentrations in the aqueous phase were assumed.To apply this assumption, for every metabolite diffusion step, the metabolite concentrations were set to the initial concentration in regions outside of the biofilm: cðx; y; z; tÞ ¼ c initial if ðx; y; zÞ is not in R biofilm Aqueous diffusion coefficients were multiplied by 0.6 for light gases and 0.25 for organic solutes to calculate the diffusion coefficient in a biofilm (Table 3) [53].Biofilm diffusion coefficients were scaled according the recommended HAL equation for PDE stable solutions [24]: where D is the metabolite diffusion coefficient (cm 2 /s), dt D is the diffusion time step (seconds), and dl is the patch length (cm) (2 μm in Table 2).Gaseous and carbon metabolite diffusion coefficients were calculated according to the respective metabolite diffusion time scales, dt D_gas and dt D_carbon , in Table 2.

Agent-based model
The agent-based model (ABM) was constructed in the established Java-based platform Hybrid Automata Library (HAL) [24].The ABM world was defined as a 3-dimensional 230 μm x 230 μm x 40 μm grid divided into 2 μm x 2 μm x 2 μm patches.Periodic boundary conditions were set in the x and y-directions.An individual agent represented an 8 μm 3 single-cell P. aeruginosa bacterium.The agent class was built using HAL's off-lattice SphericalAgent3D agent class.Each agent was assigned a unique index number.Bacteria agents were randomly initialized at z = 0 μm with a random biomass in the range of 1e -12 to 2e -12 g and a random directional angle.Biomass growth was calculated with an exponential growth rate law [5]: Where b is the updated agent biomass in grams, b 0 is the initial agent biomass in grams, μ is the growth rate informed from the agent's optimized metabolic model (hr -1 ) and dt growth is the growth time step (hr) (Table 2).If an agent biomass grew above the maximum biomass threshold (Table 4), cell division occurred, and a new daughter agent was placed in a neighboring patch of the mother agent.The mother agent biomass was randomly divided between the mother and daughter agent.The daughter agent's directional angle was set to within 10 degrees of the mother agent's angle.The daughter agent's metabolic state was set as the mother agent's metabolic state.An agent became inactive, representing cell death, in the agent-based model simulation if the predicted growth rate was 0.0 hr -1 .
Simple ABM rules for P. aeruginosa surface motility were implemented based off previous observations of physical cell-cell interactions and collective swarming movement of piliated P. aeruginosa [36,55].Pili surface motility speed was previously reported as 10 nm/s [51].Thus, for a 5 minute simulation time step, one piliated bacteria agent was estimated to move 3 μm or 1.5 patches.For each five-minute simulation time step, the surface motility ABM function was performed once for each piliated bacteria agent.To determine the agent's direction of movement, an agent on the surface (z = 0 μm) with at least one unoccupied and occupied neighbor patch set its orientation angle equivalent to a random neighbor agent.Agents moved to the new location if there were fewer than 20 cells in the surrounding region of the new location.For motility simulations, all bacteria agents in the ten-hour biofilm growth period were assumed to express pilA (S1 Fig) and perform the surface motility ABM function.For nonmotile simulations, the motility ABM function was removed for all agents.In addition, to prevent agent overlap, the built-in HAL cell mechanics algorithm was executed for each agent using force scalar and friction coefficient parameters listed in Table 3.
Each agent was assigned their own metabolic model to simulate metabolism.Four PA14 metabolic model states were available for agents to choose from, specifically, the aerobic, denitrification, denitrification +NO secretion, and oxidative stress metabolic model states.Agents decided their metabolic model state by comparing their local oxygen and NO concentrations to the respective metabolite threshold values, [O 2 ] t and [NO] t (Tables 1 and 4).For agents exposed to low oxygen and low NO environments, a stochastic parameter, R n , generated a number between 0 and 1.These cells decided between a denitrification with or without NO secretion metabolic model state by comparing their R n value to a threshold value, 0.06, which was informed from the experimental percentage of denitrification cells with a norB-limiting expression profile.After agents checked these conditions, the agents were given a metabolic state assignment in the ABM, represented by an integer value.Used to assign the correct metabolic model state to each agent in Python, this metabolic model state integer value corresponded to the order in a list of metabolic model states input into MiMICS.
A biologically-relevant 10 minute time delay, dt mRNA , for agent's to switch to a new metabolic model state upon sensing a new environmental cue was implemented, corresponding to the 10 minute delay for synthesis and detection of mRNA transcription upon a cell sensing a new environmental metabolite cue [52] (Table 2).Thus, when an agent encountered a new metabolite cue that induced a metabolic model state switch, the time the bacteria agent occupied that new metabolite environment was recorded.Once this recorded time exceeded 10 minutes [52], the bacteria agent was assigned the new metabolic model state.For dividing agents, this recorded time was equally split between the mother and daughter agents, representing observations of detected mRNA amounts in dividing cells [52].

Coupled model components in MiMICS
For each simulation time step, dt sim , MiMICS executed each model component in the following order: 1) perform agent functions in the ABM, 2) constrain and optimize agent metabolic models, and 3) solve metabolite reaction-diffusion PDEs.A bacterial growth time scale was estimated as 30 minutes based on the doubling time of P. aeruginosa in SCFM [29].Bacterial growth was determined as the rate-limiting biological process.Due to the expected stochasticity of bacteria agent division events, one MiMICS simulation time step, dt sim , was set as 5 minutes.
At the beginning of the time step, agent functions were performed in the ABM in order of agent death, division, metabolic state assignment, motility, and agent mechanics.Next, each agent's biomass, index value, metabolic state value, and the agent's local metabolite concentrations (i.e.oxygen, nitrate, NO, glucose) was passed from the ABM (Java) to the GENRE (Python).The metabolite concentrations were converted to fluxes using Eq 1 to constrain the agent's metabolic model exchange fluxes.In transcriptome-guided MiMICS simulations, the agent's metabolic state assignment number was used for agents to select a metabolic model state from the four PA14 metabolic model states.In transcriptome-free MiMCIS simulations, all agents simulated metabolism with PA14 GENRE constrained on SCFM, but not integrated with transcriptomic data.Each agent's GENRE was optimized for a biomass growth rate that was used to calculate (Eq 5) and update the agent biomass in the ABM.In addition, each agent's optimized metabolic model predicted metabolite secretion and consumption fluxes.The oxygen, nitrate, NO, glucose exchange fluxes predicted by each agent's metabolic model were converted to concentrations (Eq 1) and passed to the metabolite PDEs to update the metabolite concentrations at the metabolite grid location where the agent resided using Eq 2.Last in the simulation time step, metabolite diffusion was simulated by solving the continuum-scale PDE models for each metabolite grid (i.e.oxygen, nitrate, NO, glucose).
Metabolite consumption and production time scales were estimated as 0.05 seconds based on metabolite maximum uptake kinetics of oxygen [2].A general metabolite diffusion time scale, dt D , was estimated as 0.02 seconds based on metabolite diffusion across a one-dimensional patch with length 2 μm.Compared to the biomass growth time scale, metabolite reaction and diffusion processes were considered to be at pseudo-steady state, in agreement with previous multi-scale models [59].
When metabolite PDEs were solved, metabolite diffusion was simulated for multiple time steps for every metabolic reaction time step.This approach reconciled the difference between the relatively slow metabolic reaction time scale (0.05s) compared to the fast metabolite diffusion time scale (0.02s), while maintaining stable PDE solutions.The number of times diffusion was performed for carbon and gaseous metabolites was calculated by dt ex /dt D_carbon and dt ex / dt D_gas , respectively.Accordingly, for carbon substrates with lower magnitude diffusion coefficients and dt D carbon similar in magnitude to dt ex , metabolite diffusion was performed 5 times per five-minute simulation time step (Table 2).To achieve stable diffusion coefficients for gases with relatively high magnitude diffusion coefficients and dt D gas lower in magnitude compared to dt ex , metabolite diffusion was performed 20 times per five-minute simulation time step.
An experimentally observed growth lag phase was estimated as ~1.5 hours [5] which translated to ~20 five-minute simulation time steps.Therefore, MiMICS simulations were run for 100 five-minute simulation time steps, which represented ten hours of PA14 biofilm growth accounting for the experimentally observed growth lag phase (Table 2).At desired time points, MiMICS output each agent's index, biomass, growth rate, metabolic state assignment number, intracellular metabolic reaction fluxes, spatial coordinates, and local metabolite concentrations.In addition, the metabolite concentrations for each metabolite grid were output.

Parameterization
The oxygen threshold parameter, [O 2 ] t , was varied from 0.19-0.23 mM.The mean absolute error between experiment and MiMICS outputs was calculated from 21 replicate simulations performed for each [O 2 ] t parameter value.The total MiMICS model error was calculated from the summation of the error from MiMICS predictions between two experimental outputs: the percentage of denitrification classified cells, and the percentage of oxidative stress classified cells.Denitrification cells were classified to have expression of one or more the following denitrification genes: nirS, norB, nosZ, and no expression of the oxidative stress katA gene.Oxidative stress cells were classified to have expression of the oxidative stress katA gene.The [O 2 ] t parameter value with the smallest error relative to the experiment was used for subsequent simulations.

MiMICS specifications and code availability
In MiMICS, a Py4J JavaGateway Server was implemented to interface the Java-based HAL platform and Python-based metabolic model.Through this Py4J JavaGateway, functions are called from Python to initialize and run the ABM.Additionally, Python functions called through the JavaGateway retrieve and pass information between Python and Java and save spatial agent and metabolite information at each simulation time point.MiMICS was initially developed on a Mac computer (macOS Catalina v.10.15.7, 1.1 GHz Dual-Core Intel Core i3) using IntelliJ IDEA and 1-2 central processing units (CPUs).All MiMICS simulations were executed on Rivanna High-Performance Computing (HPC) System at the University of Virginia.MiMICS Java files were compiled as a JAR file to run on the Rivanna HPC System.To reduce simulation runtime, MiMICS split metabolic model calculations for each agent across 35 CPUs using Python multiprocessing v.0.70.14.One MiMICS simulation of the ten-hour biofilm growth period on Rivanna HPC System across 35 CPUs took 40 minutes and required 40 GB of memory.When appropriate, multiple MiMICS simulations (e.g. for simulation replicates, parameter sweeps) were executed across multiple computing nodes to prevent multiple JavaGateway servers running on a single computing node.MiMICS simulations used Python v.3.9, cobra v.0.29.0, py4j v.0.10.9.7, Java v.1.8.0, HAL v.1.1.0.MiMICS source code and detailed user guide can be found at: https://github.com/tracykuper/mimics.

Image analysis and data visualization
Biofilm image analysis of MiMICS simulation outputs with and without motility was performed using Python's scikit-image (v.0.19.0) regionprops function to measure total biofilm area and average cluster area.Experimental and simulation data and 2D images were plotted with the matplotlib python package (v.3.6.3).3D images were plotted with plotly package (v.4.14.3)

Statistical analysis
Statistical analysis was performed with the scipy stats python package (v.1.10.0).Variation among biofilm structure metrics among experimental and simulation conditions was performed with analysis of variance.Pairwise comparisons of individual group means were performed using a Tukey post hoc analysis.Values of p < 0.05 were considered statistically significant.

Spatial neighborhood analysis
Bulk neighborhood analysis of simulation data was performed similarly to the experimental neighborhood analysis [5].For each gene of interest, agents with active flux encoded by the reaction in the 99 th percentile were selected as "center" agents.Using the 3D centroid agent coordinates, distances between center agents and nearby agents within 3 μm were calculated.The closest five agents to the center agent were selected as neighbor agents.Neighbor agents were grouped together and the average reaction flux of the neighbors was calculated.This average neighbor reaction flux was divided by the corresponding population average reaction flux, omitting the center agents in the population.For each pair of genes, a Pearson correlation was calculated to determine spatially correlated genes among 50 simulation replicates.

Fig 1 .
Fig 1. Overview of the MiMICS computational framework.Individual bacterial agents were initialized at t = 0hrs in a 3D ABM world.Metabolites were initialized based on user-defined metabolite concentrations.For each five-minute simulation time step, the MiMICS framework was executed, consisting of a set of-omics data-integrated metabolic models, an ABM, and a metabolite reaction-diffusion model.MiMICS simulation outputs included extracellular metabolite concentrations, as well as each agent's location, metabolic state, and intracellular metabolite fluxes at each time point.https://doi.org/10.1371/journal.pcbi.1012031.g001

Fig 2 .
Fig 2. MiMICS incorporated four unique spatial-transcriptome-guided PA14 metabolic model states.(A) Heatmap of experimental gene expression values of each PA14 metabolic model state.(B) Heatmap of predicted flux values of gene-encoded reactions for each PA14 metabolic model state.Flux values reported are absolute flux values normalized to the maximum absolute flux value among metabolic model states.(C) PA14 GENRE denitrification and oxidative stress pathway updated from iPau21 GENRE with NO secretion pathway, NO iron cytotoxic reaction, and katA-encoded oxidative stress reaction.(D) Predicted fluxes of metabolites exchanged with the extracellular environment for a transcriptome-free PA14 metabolic model and the four transcriptome-guided PA14 metabolic model states.Fluxes reported are normalized to the maximum flux among metabolic models.https://doi.org/10.1371/journal.pcbi.1012031.g002 Fig 3).Consequently, reflective of the intracellular metabolic model states (Fig 2B), transcriptome-guided MiMICS predicted a relatively low number of cells (10% of cells) utilizing narG-and napAencoded reactions compared to the cell population utilizing norB-and nosZ-encoded reactions (15% of cells) (Fig 4).This MiMICS prediction agreed with the experimental results reporting that relatively few cells expressed narG and napA, 1.5% and 4.3% of cells, respectively, compared to cell populations with high norB and nosZ expression, 12.1% and 11.8% of cells, respectively (Fig 4).

Fig 3 .
Fig 3. Transcriptome-guided MiMICS improved predictions of microscale metabolic heterogeneity and NO microenvironment in PA14 biofilm.Shown are representative 3D renderings of PA14 biofilms grown for ten-hours from the experiment, and transcriptome-free and transcriptome-guided MiMICS simulations.Cells are colored according to their metabolic state.Shown are 2D yz slices (x = 160 μm, 22 μm, 100 μm for experimental, transcriptome-guided MiMICS, and transcriptome-free MiMICS, respectively) of cell metabolic states, and predicted oxygen and NO concentrations.The x-values were chosen to compare similar biofilm colony structures across experimental and simulation conditions.Experimental data was reconstructed from Dar and co-workers.https://doi.org/10.1371/journal.pcbi.1012031.g003

Fig 4 .
Fig 4. Transcriptome-guided MiMICS improved prediction accuracy of heterogeneous intracellular metabolism in PA14 biofilm compared to transcriptome-free MiMICS.From experimental data, plotted is the average percentage of total cells expressing a gene.From transcriptome-guided MiMICS and transcriptome-free MiMICS simulation data, plotted are the average percentage of total agents with active flux through the reaction encoded by a gene.Error bars represent one standard deviation from n = 7 experimental replicates and n = 50 simulation replicates.Experimental data provided by Dar and coworkers.https://doi.org/10.1371/journal.pcbi.1012031.g004 (Fig 5B).Genes were compared which had both positive experimental expression values and simulated active flux values of reactions encoded by the respective gene.Shown in Fig 5B, MiMICS accurately predicted the spatial correlation for 33 gene pairs, and incorrectly predicted the spatial correlation for 22

Fig 5 .Fig 6 .
Fig 5. MiMICS captured spatial relationships of intracellular metabolism in PA14 biofilm.(A) Representative xy projections of PA14 biofilms from experiments and transcriptome-guided MiMICS simulations.Cells plotted are located near the z = 0 μm surface.Colored cells have high expression of the gene listed (experiment) or high reaction flux encoded by the gene listed (simulation).Circled areas highlight regions of interest where there is an anticorrelation between TCA cycle metabolism with denitrification and oxidative stress metabolism.Scale bar represents 20 μm.(B) Neighborhood gene spatial correlation comparison between experiment and transcriptome-guided MiMICS simulation.Spatial correlation between gene pairs was assessed by a Pearson correlation, where +1 and -1 value correspond to a strong positive and strong negative spatial correlation, respectively.The experimental and simulation Pearson correlation values are plotted in the upper right and lower left triangle of each square, respectively.Simulation Pearson correlation values were determined from 50 simulation replicates.Experimental data was reconstructed from Dar and co-workers.https://doi.org/10.1371/journal.pcbi.1012031.g005